\documentclass[twoside,a4paper]{article}
\usepackage{geometry}
\geometry{margin=1.5cm, vmargin={0pt,1cm}}
\setlength{\topmargin}{-1cm}
\setlength{\paperheight}{29.7cm}
\setlength{\textheight}{25.3cm}

% useful packages.
\usepackage[UTF8]{ctex}	% 中文支持
\usepackage{listings} 	% 插入代码
\usepackage{xcolor} 	% 语法高亮
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{enumerate}
\usepackage{graphicx}
\usepackage{multicol}
\usepackage{fancyhdr}
\usepackage{layout}

% some common command
\newcommand{\dif}{\mathrm{d}}
\newcommand{\avg}[1]{\left\langle #1 \right\rangle}
\newcommand{\difFrac}[2]{\frac{\dif #1}{\dif #2}}
\newcommand{\pdfFrac}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\OFL}{\mathrm{OFL}}
\newcommand{\UFL}{\mathrm{UFL}}
\newcommand{\fl}{\mathrm{fl}}
\newcommand{\op}{\odot}
\newcommand{\Eabs}{E_{\mathrm{abs}}}
\newcommand{\Erel}{E_{\mathrm{rel}}}

\newcommand{\R}{\mathbb{R}}
\renewcommand{\proofname}{Proof.}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}{Lemma}[section]
\newtheorem{exercise}{Exercise}[section]
\newtheorem{definition}{Definition}[section]

\begin{document}
	
	\pagestyle{fancy}
	\fancyhead{}
	\lhead{行一凡 (3190105815)}
	\chead{Final Project}
	\rhead{\today}
	
	\section{项目结构}
	程序所需头文件在 include 文件夹中，主程序文件在主目录下；三个文件夹 matlab、outputs、pictures 依次存放 matlab 绘图脚本、程序输出结果、生成图片。进入主目录下
	\begin{itemize}
		\item 输入 make story 编译生成 report.pdf 文件；
		\item 输入 make test 命令编译所有测试程序；
		\item 输入 make run 运行程序，将结果输出到指定文件当中；
		\item 输入 make clean 清除可执行程序。
	\end{itemize}
	
	\section{程序结构}
	\begin{enumerate}[1.]
		\item include 文件中包含了程序所需头文件的相关代码。
		\begin{itemize}
			\item Tensor.h 主要定义了 Tensor 类，作为基本储存容器；同时定义了所需的宏、变量别名和常用的简易函数。Func、FuncX 和 FuncXt 定义了一元向量函数、多元向量函数和将时间 t 作为单独参数的多元向量函数；
			
			\item Polynomial.h 定义了多项式类，用于计算插值多项式；
			
			\item Integrator.h 定义了区域 Domain 和积分算子 Integrator，分别用于标记积分域和计算体积分、面积分；
			
			\item GhostFiller.h 定义了填充器，对于指定指标位置获取周围的单元格点或 ghost cell；
			
			\item FVOperator.h 定义了网格插值和约束算子和有限体积法所需的 Laplacian、Gradient 和 Divergence 算子；
			
			\item RungeKutta.h 定义了隐式-显式 RK 方法类，用于存放迭代系数和权；
			
			\item MultigridSolver.h 定义了多重网格法，使用 V-循环和 G-S 迭代法求解；
			
			\item Projection.h 定义了投影算子 $ P=I-GL^{-1}D $；
			
			\item AdvectionDiffusion.h 定义了对流扩散方程的求解器；
			
			\item INSE.h 定义了 INSE 的求解器。
		\end{itemize}
		
		\item 主目录中包含测试文件
		\begin{itemize}
			\item test1.cpp 测试 Traveling sinusoidal waves；
			
			\item test2.cpp 测试 Gaussian patch in vortex shear；
			
			\item test31.cpp 和 test32.cpp 测试 Taylor vortex。
		\end{itemize}
	\end{enumerate}
	 
	
	\section{数值实现}
	\subsection{数值积分}
	先考虑一维积分。取
	\begin{equation}\label{key}
		f_i(x) = \begin{cases}
			1, & x = i\\
			0, & x\neq i
		\end{cases},\quad i = 0,1,\cdots,n
	\end{equation}
	以及 $ x_i=0,1,\cdots,n $，其中 $ n+1 $ 是取点的个数。通过计算 $ f_i $ 在 $ x_i $ 上的插值多项式积分计算格点函数值的权
	\begin{equation}\label{key}
		w_i = \int_0^n p_i(x)dx
	\end{equation}
	于是对于给定区间 $ [a,b] $ 和被积函数 $ f $，如果均匀取 $ n+1 $ 个点，则
	\begin{equation}\label{key}
		\int_a^bf(x)dx \approx h\sum_{i=0}^{n}w_if(x_i),\quad x_i = a + hi,\ h=\frac{b-a}{n}
	\end{equation}
	对于 $ Dim $ 维的情形，利用多重积分，在每个方向上叠加权重。例如 $ Dim=2 $ 时，权系数为 $ W_{ij} = w_iw_j $，即 $ W $ 为 2 阶张量。对于给定方形区域 $ S =[a,b]\times [c,d] $ 和被积函数 $ f(x,y) $，有
	\begin{equation}\label{key}
		\iint_S f(x,y)dxdy \approx h^2\sum_{i=0}^{n}\sum_{j=0}^{n}W_{ij}f(x_i,y_j),\quad x_i = a + hi,\ y_j=c_i+hj
	\end{equation}
	更高维的情形同理提高张量阶数即可。
	
	\subsection{单元填充器}
	有限体积的算子作用需要周围的格点，这时就需要利用填充器获取当前单元周围存在或假定存在的格点的值。通常需要每个方向左右各 2 个格点，为了便于索引，创建并返回一个长度为 $ 5 $ 的容器。
	
	对于周期边界，在 Tensor 类中定义特定的函数，使得我们可以抵达一个边界时将索引返回到对面的边界，从而简化处理；对于 Dirichlet 和 Neumann 边界，通过 GhostFiller 携带的边界函数计算并填充。
	
	对于非周期边界条件下的特定的算子，例如 Laplacian，在求解过程中 $ I-k\nu\gamma L $ 作为系数矩阵，需要将边界条件分离。例如
	\begin{equation}\label{key}
		\left<\phi\right>_{\mathbf{i}+\mathbf{e}^d} = \frac{1}{3}(-13\left<\phi\right>_{\mathbf{i}}+5\left<\phi\right>_{\mathbf{i}-\mathbf{e}^d}-\left<\phi\right>_{\mathbf{i}-2\mathbf{e}^d})
	\end{equation}
	其中 Dirichlet 边界条件 $ 4\left<g\right>_{\mathbf{i}+\frac{1}{2}\mathbf{e}^d} $ 无法作为系数矩阵的一部分，因此在计算 Laplacian 的 ghost cell 时将其排除。缺少的边界条件通过 GhostFiller 补充到线性系统的右端。
	
	\subsection{有限体积算子}
	网格插值和约束算子和 Laplacian、Gradient 的计算方法利用填充器和 4 阶近似公式计算。对于 Divergence 算子，在计算 $ \left<u_d\phi\right>_{\mathbf{i}} $ 时应当区分周期边界和非周期边界。周期边界条件下应当利用
	\begin{equation}\label{key}
		\begin{aligned}
			\mathbf{F}\left<u_d,\phi\right>_{\mathbf{i}+\frac{1}{2}\mathbf{e}^d} &= \left<u_d\phi\right>_{\mathbf{i}+\frac{1}{2}\mathbf{e}^d}\\
			\mathbf{D}\left<\mathbf{u}\phi\right>_{\mathbf{i}} &= \frac{1}{h}\sum_{d=1}^{D}(\left<u_d\phi\right>_{\mathbf{i}+\frac{1}{2}\mathbf{e}^d}-\left<u_d\phi\right>_{\mathbf{i}-\frac{1}{2}\mathbf{e}^d})
		\end{aligned}
	\end{equation}
	直接依次代入计算；而对于非周期边界，在取 ghost cell 时需要分别计算 $ u_d $ 和 $ \phi $ 的边界条件。
	
	投影算子 $ P=I-GL^{-1}D $ 按照其中每部分算子顺序作用。对于 $ L^{-1} $，在周期边界上求解 Poisson 方程即可。
	
	\subsection{隐式-显式 Runge-Kutta 方法}
	使用隐式-显式 Runge-Kutta 方法求解
	\begin{equation}\label{key}
		\pdfFrac{\phi}{t} = -\nabla\cdot (\mathbf{u}\phi) + \nu\Delta\phi + f
	\end{equation}
	利用有限体积法，求解
	\begin{equation}\label{key}
		\begin{aligned}
			\difFrac{\phi}{t} &= \mathbf{X}^{[\mathrm{E}]}(\phi,t) + \mathbf{X}^{[\mathrm{I}]}(\phi)\\ \mathbf{X}^{[\mathrm{E}]}(\phi,t) &= \left<f\right> - \mathbf{D}\left<\mathbf{u}\phi\right>\\
			\mathbf{X}^{[\mathrm{I}]}(\phi) &= \nu \mathbf{L}\left<\phi\right>
		\end{aligned}
	\end{equation}
	在每一步使用多重网格法求解
	\begin{equation}\label{key}
		(I-k\nu\gamma L)\left<\phi\right>^{(i)} = \left<\phi\right>^n + k\sum_{j=1}^{i-1}a_{i,j}^{[\mathrm{E}]}\left[\left<f\right>^{(j)} - D\left<\mathbf{u}\phi\right>^{(j)}\right] + k\nu\sum_{j=1}^{i-1}a_{i,j}^{[\mathrm{I}]}L\left<\phi\right>^{(j)}
	\end{equation}
	由于在 $ I-k\nu\gamma L $ 作用时分离了边界条件，因此需要在右端补充边界条件。最后加权推进
	\begin{equation}\label{key}
		\left<\phi\right>^{n+1} = \left<\phi\right>^n + k\sum_{j=1}^{s-1}b_{s}^{[\mathrm{E}]}\left[\left<f\right>^{(j)} - D\left<\mathbf{u}\phi\right>^{(j)}\right] + k\nu\sum_{j=1}^{s-1}b_{s}^{[\mathrm{I}]}L\left<\phi\right>^{(j)}
	\end{equation}
	
	求解不可压 NS 方程的方法只需要在上述过程中添加投影算子 $ P $ 的作用即可。最后在特定时间，利用
	\begin{equation}\label{key}
		\mathbf{L}\left<p\right> = \mathbf{D}(\left<\mathbf{g}\right>-\mathbf{D}\left<\mathbf{uu}\right>+\nu \mathbf{L}\left<\mathbf{g}\right>)
	\end{equation}
	在周期边界上求解 Poisson 方程可得到 $ p $。由于周期边界条件下 Laplacian 算子奇异，因此得到解 $ \left<p\right> $ 后还需要减去与真实解的差得到正确的计算解。
	
	\section{测试内容}
	\subsection{Traveling sinusoidal waves}
	对于 $ \nu=0.01,\ \mathbf{u} = (1.0,0.5),\ \mathbf{k}=(2\pi,4\pi),\ \theta = (k_1x-u_1t,k_2x-u_2t) $，精确解 $ \phi(x,y,t) = \sin\theta_1\sin\theta_2 $，压力项为
	\begin{equation}\label{key}
		f(x,y,t) = \nu(k_1^2+k_2^2)\sin\theta_1\sin\theta_2 + u_1(k_1-1)\cos\theta_1\sin\theta_2 + u_2(k_2-1)\cos\theta_2\sin\theta_1
	\end{equation}
	测试 $ h = \frac{1}{16},\frac{1}{32},\frac{1}{64},\ \mathrm{Cr} = 1 $，在 $ [0,1]^2 $ 上以下边界为 Dirichlet 条件，上边界为 Neumann 条件，从 $ t_0 $ 推进到 $ t_e=1 $。具体计算数据如下：
	\begin{center}
		\begin{tabular}{|c|c|c|}
			\hline
			N & 误差范数 & 运行时间(s)\\
			\hline
			16 & 0.0144725 & 77.072\\
			32 & 0.00110599 & 396.022\\
			64 & 7.84959e-005 & 2955.74\\
			\hline
		\end{tabular}
	\end{center}
	根据表格数据，可以得出求解算法平均收敛阶数为 $ 3.7632 $，符合预期的收敛阶数。 
	
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.55\linewidth]{pictures/test1-16}
		\label{fig:test1-16}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.55\linewidth]{pictures/test1-32}
		\label{fig:test1-32}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.55\linewidth]{pictures/test1-64}
		\label{fig:test1-64}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.55\linewidth]{pictures/test1-convergence}
		\caption{误差范数取 $ 2 $ 的对数}
		\label{fig:test1-convergence}
	\end{figure}

	\newpage
	\subsection{Gaussian patch in vortex shear}
	此测试的速度场为
	\begin{equation}\label{key}
		\mathbf{u}(x,y) = a_V(\sin^2(\pi x)\sin(2\pi y),-\sin(2\pi x)\sin^2(\pi y))
	\end{equation}
	其中 $ a_V = 0.1 $。我们使用精确解 $ \phi(x,y,t) = \sin\theta_1\sin\theta_2 $ 来计算初值，$ \theta=(2\pi x-t, 4\pi y-t) $，推导得到压力项
	\begin{equation}\label{key}
		\begin{aligned}
			f(x,y,t) &= 20\nu\pi^2\sin\theta_1\sin\theta_2 -\cos\theta_1\sin\theta_2 -\cos\theta_2\sin\theta_1\\
			&+2a_V\pi\phi(x,y,t)[\sin(\pi x)\cos(\pi x)\sin(2\pi y)-\sin(2\pi x)\sin(\pi y)\cos(\pi y)]\\
			&+ a_V[2\pi\sin^2(\pi x)\sin(2\pi y)\cos\theta_1\sin\theta_2 - 4\pi\sin(2\pi x)\sin^2(\pi y)\sin\theta_1\cos\theta_2]
		\end{aligned}
	\end{equation}
	测试 $ \nu=0.001,\ h = \frac{1}{16},\frac{1}{32},\frac{1}{64},\ \mathrm{Cr} = 1 $，在 $ [0,1]^2 $ 上应用周期边界条件，从 $ t_0 $ 推进到 $ t_e=1 $。精确解通过 Richardson 外推法计算。具体计算数据如下：
	\begin{center}
		\begin{tabular}{|c|c|c|}
			\hline
			N & 误差范数 & 运行时间(s)\\
			\hline
			16 & 1.52313e-007 & 39.358\\
			32 & 1.04233e-008 & 244.565\\
			64 & 6.58599e-010 & 2283.88\\
			\hline
		\end{tabular}
	\end{center}
	其中运行时间不包括外推法计算精确解的时间。根据表格数据，可以得出求解算法平均收敛阶数为 $ 3.9267 $，符合预期的收敛阶数。 
	
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.5\linewidth]{pictures/test2-16}
		\label{fig:test2-16}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.5\linewidth]{pictures/test2-32}
		\label{fig:test2-32}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.55\linewidth]{pictures/test2-64}
		\label{fig:test2-64}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.55\linewidth]{pictures/test2-convergence}
		\caption{误差范数取 $ 2 $ 的对数}
		\label{fig:test2-convergence}
	\end{figure}
	
	\subsection{Taylor vortex}
	使用精确解
	\begin{equation}\label{key}
		\begin{aligned}
			\mathbf{u}(x,y,t) &= 1+2\exp(-8\pi^2\nu t)
			\begin{pmatrix}
				-\cos(2\pi(x-t))\sin(2\pi(y-t))\\
				\sin(2\pi(x-t))\cos(2\pi(y-t))
			\end{pmatrix}\\
			p(x,y,t) &= -\exp(-16\pi^2\nu t)(\cos(4\pi(x-t))+\cos(4\pi(y-t)))
		\end{aligned}
	\end{equation}
	测试 $ \mathrm{Cr} =0.75,\ \mathrm{Re}=30,300,\ h = \frac{1}{16},\frac{1}{32},\frac{1}{64} $，在 $ [0,1]^2 $ 上应用周期边界条件，从 $ t_0 $ 推进到 $ t_e=0.5 $。
	
	\begin{table}[!htb]
		\centering
		\begin{tabular}{|c|c|c|c|}
			\hline
			N & 误差范数 $ (u) $ & 误差范数 $ (p) $ & 运行时间(s)\\
			\hline
			16 & 0.0149588 & 0.0104219 & 68.664\\
			32 & 0.000531599 & 0.000538536 & 532.13\\
			64 & 1.68808e-005 & 2.48418e-005 & 6652.34\\
			\hline
		\end{tabular}
		\caption{$ \mathrm{Re} = 30 $}
	\end{table}
		
	\begin{table}
		\centering
		\begin{tabular}{|c|c|c|c|}
			\hline
			N & 误差范数 $ (u) $ & 误差范数 $ (p) $ & 运行时间(s)\\
			\hline
			16 & 0.301677 & 0.414178 & 68.901\\
			32 & 0.00840476 & 0.0195737 & 532.739\\
			64 & 0.000299365 & 0.000758481 & 6344.96\\
			\hline
		\end{tabular}
		\caption{$ \mathrm{Re} = 300 $}
	\end{table}
	根据表格数据，可以得出 $ \mathrm{Re}=30,300 $ 两种情况下，对速度场 $ u $ 的平均收敛阶数为 $ 4.8957 $ 和 $ 3.4577 $，对压力项 $ p $ 的平均收敛阶数为 $ 4.3563 $ 和 $ 3.1514 $，基本符合预期的收敛阶数。 
	
	由于两种情况下的解类似，只展示 $ \mathrm{Re}=30 $ 的情况。下图将速度 $ u $ 分解为 $ x,y $ 两部分显示：
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.6\linewidth]{pictures/test31-16x}
		\label{fig:test31-16x}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.6\linewidth]{pictures/test31-16y}
		\label{fig:test31-16y}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.55\linewidth]{pictures/test31-16p}
		\label{fig:test31-16p}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.55\linewidth]{pictures/test31-32x}
		\label{fig:test31-32x}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.55\linewidth]{pictures/test31-32y}
		\label{fig:test31-32y}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.55\linewidth]{pictures/test31-32p}
		\label{fig:test31-32p}
	\end{figure}
	
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.55\linewidth]{pictures/test31-64x}
		\label{fig:test31-64x}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.55\linewidth]{pictures/test31-64y}
		\label{fig:test31-64y}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.55\linewidth]{pictures/test31-64p}
		\label{fig:test31-64p}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.55\linewidth]{pictures/test31-convergence}
		\caption{误差范数取 $ 2 $ 的对数}
		\label{fig:test31-convergence}
	\end{figure}
	
	
	

\end{document}